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Abstract. - We report on the nature of the thermal denaturation transition of homogeneous 
DNA as determined from a renormalisation group analysis of the Peyrard-Bishop-Dauxois model. 
Our approach is based on an analogy with the phenomenon of critical wetting that goes further 
than previous qualitative comparisons, and shows that the transition is continuous for the average 
base-pair separation. However, since the range of universal critical behaviour appears to be very 
narrow, numerically observed denaturation transitions may look first-order, as it has been reported 
in the literature. 
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Introduction. — In addition to its central relevance 
in biology, the DNA molecule also displays a variety of re- 
markable physical properties, some of which are key to the 
understanding of DNA function [I]. For instance, mechan- 
ical properties such as bending, twisting, or compression 
are directly related to DNA replication or transcription, 
which requires that the two strands are separated in or- 
der that they can be read by DNA or RNA polymerase. 
This can be achieved by various mechanisms, including 
pulling enzymes, mechanical force or gentle heating. In 
this latter case, the process is known as DNA thermal 
denaturation or DNA melting, and has received a great 
deal of attention over several decades [2]. Experimental 
observation of the fraction of bound pairs or the average 
base-pair separation as a function of temperature reveals 
a sharp jump in the denaturation curves from double- to 
single-stranded DNA, a behaviour that hints at some sort 
of phase transition. However, controversies remain regard- 
ing the nature of this transition (whether first or second 
order), see below. This point goes beyond mere physical 
curiosity and has biological relevance, for there is increas- 
ing interest in the correspondence between functional and 
thermodynamic melting properties, e.g., the identification 
of coding sequences in genomes on the basis of thermody- 



namic melting behaviour [3]. 

Models of varying complexity and applicability have 
been developed that account for DNA melting. Two 
large families are based on the models of Poland-Scheraga 
(PS) and Peyrard-Bishop-Dauxois (PBD) [5^, or mod- 
ifications thereof. Within the PS framework, the DNA is 
described as a sequence of base pairs that can be either 
bound or unbound. Thermal fluctuations cause segments 
of DNA to unbind, creating temporarily denaturated loops 
of variable size which can ultimately coalesce upon in- 
creasing the temperature, thus triggering the denatura- 
tion transition. It has been shown that the entropy con- 
tribution of the loops depends on their size, I, as ^ l/^*^, 
and three different scenarios have been reported depend- 
ing only on the value of c: c < 1, no phase transition, 
1 < c < 2, continuous transition, and c > 2, first-order 
phase transition. The value of c is not easily determined, 
but it has been recently demonstrated that taking into 
account the excluded-volume interactions between denat- 
urated loops and the rest of the chain is enough to give 
c > 2 [6^, so the transition is therefore first-order (see [7] 
for a refined analysis of these issues and a discussion of 
some open questions). 

Turning to PBD-type models, the situation is much less 
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clear. The PBD model considers only the stretching be- 
tween corresponding base-pairs (see details next section). 
The transition proceeds as described for the PS, but in- 
cluding intermediate states because the stretching is a con- 
tinuously varying variable. It has been thoroughly studied 
by means of Monte Carlo simulations, Langevin dynam- 
ics, path integral methods and different transfer integral 
approaches, but the question of the order of the transition 
remains as yet unsettled. Claims have been reported in 
the literature that the transition is first-order yet with a 
diverging correlation length, asymptotically second-order 
although very sharp looking in appearance, while other 
studies are inconclusive. We shall return to this point in 
the next section. 

In this Letter we take on the question of the order of 
the DNA denaturation transition in the PBD model by 
means of an exact renormalisation group analysis. Our 
approach is based on an analogy with the phenomenon of 
critical wetting that goes further than previous qualitative 
comparisons. In the next section, we review the PB and 
PBD models. Next, we explain the analogy with wetting 
and then proceed with the renormalisation group calcu- 
lation, which shows that the transition is continuous for 
the average base-pair separation. A summary with our 
conclusions is presented in the final section. 

The Peyrard-Bishop-Dauxois model of DNA. 

The Peyrard-Bishop model ignores the helicoidal structure 
of the DNA molecule and the properties associated with 
it, and focuses on the stretching of the hydrogen bonds 
connecting base pairs, which are represented by continu- 
ous variables hn {n = 1,2, N, where N is the length 
of the chain). For homogeneous samples (only AT of GC 
pairs), the Hamiltonian of the model reads [5] 



N 



+ W{K,h,,^i) + V{K) 



(1) 



where the first term is the kinetic energy for bases of 
mass m, W{hn, hn-i) — k{hn — /in-i)^ describes the har- 
monic stacking interaction between neighbouring bases, 
and V represents the average potential between the two 
bases in a pair which is modeled by a Morse potential, 
V{hn) = D{e~^^'^ — 1)^. Note that the asymmetry of the 
two strands is neglected in that a common mass m for the 
bases is used and the same stacking coefficient k along the 
chain is assumed. D is the dissociation energy of the pairs 
and a denotes the spatial range of the potential. Their 
precise values, which are unimportant for our purposes, 
can be determined from the fitting of DNA experimental 
denaturation curves ^ . Two standard observables are the 
average stretching {h) and the density of bound base pairs 

It was soon recognised that the simple PB model needed 
to be improved if it was to properly account for the sharp 
shape of the experimental denaturation curves. This was 
achieved by changing from a harmonic stacking interaction 



to the nonharmonic one 



(2) 



whose origin lies in the change in the electronic distribu- 
tion on the bases when the hydrogen bonds are broken and 
that provides a more realistic treatment of the phosphate 
backbone stiffness. This new term leads to very sharp 
melting transitions at substantially reduced denaturation 
temperatures [S]. 

Let us now briefiy review the different scenarios that 
have been reported for the denaturation transition de- 
pending on the stiffness parameter p, for both homo- 
geneous and heterogeneous DNA samples. There is a 
consensus that, in the simplest case p — Q, second- 
order denaturation transitions are observed irrespective 
of the composition of the sample (whether heterogeneous 
or homogeneous) |10[|11| . A nonvanishing p and het- 
erogeneous sequences successfully exhibit the character- 
istic abrupt, multistep melting observed in heterogeneous 
DNA molecules pjl for intermediate-length sequences (see 
also [H]). On the contrary, in the case of nonzero p and 
homogeneous DNA the situation is much less clear. On the 
one hand, there is a claim by Cule et al. that the transi- 
tion region is extremely narrow, making it very sharp in 
appearance although, asymptotically, it is expected to be 
second-order [TU]. On the other hand, Dauxois et al. [5] re- 
ported a first-order transition, yet with a diverging corre- 
lation length [111. Subsequent improved transfer-integral 
investigations by Joyeux et al. [13] did not settle the ques- 
tion, inasmuch as the numerics seems to indicate a first- 
order transition, but without discarding the possibility of 
a narrow second-order one. Finally, results from a very 
recent path integral investigation suggest that the denat- 
uration of homogeneous DNA has the features of a second- 
order phase transition |14j . In what follows, we clarify this 
situation by exploiting an analogy with the phenomenon 
of wetting. 

The wetting analogy. — The canonical partition 
function of the model factorises as usual into a prod- 
uct of kinetic and configurational parts, Z = ZpZy, with 
Zp = [2TTmkBT)^/'^ and 



N 



W dh„ 



(3) 



n=0 



where we have defined H' ~ {W + V) as the configu- 
rational part of H . Given that most experiments on DNA 
thermal denaturation are performed in water, the kinetic 
term does not play any role, and hence we can restrict 
ourselves on the configurational part. In the continuum 
limit, for small values of hn — hn-i, the configurational 
part of the PBD Hamiltonian can be expressed as 



dx 



(1 -I- pe-^°'"){Vhy + wie-"" + wae 



-2ah 



(4) 
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where wi, 1x12 and k are generic parameters. For p = 
TJetu is the standard interfacial Hamiltonian for equihb- 
rium critical wetting transitions in the presence of short- 
ranged forces, that is, the unbinding of the interface sep- 
arating two coexisting phases from a substrate, which oc- 
curs upon increasing the temperature |15| . It constitutes 
an approximation to the PBD Hamiltonian that disregards 
the kinetic terms, but from which equilibrium information 
can be gleaned. The two strands of the DNA molecule 
correspond to the substrate and the interface in the wet- 
ting context, and the denaturation of the former to the 
unbinding of the latter. The analogy also extends to het- 
erogeneous sequences, the thermal denaturation of het- 
erogeneous DNA corresponding to the wetting of a one- 
dimensional interface from a disordered substrate. This 
formal relation between wetting and DNA thermal denat- 
uration was already noticed by Fisher [TB] and exploited by 
Cule and Hwa [TU] and by Ares et al. [T7] . The analogy was 
carried further in [18], where it was pointed out that the 
reported critical exponents characterising the DNA denat- 
uration transition in the homogeneous case, {h) ~ |<5|'' and 
£, ~ \5\-'' (where 5 = [T - Tc)/Tc and ^ is the correlation 
length) [TT] , are those of two-dimensional critical wetting, 
/3 = — 1 and v = 2 [15j . It was also shown by numer- 
ical simulations that the average stretching {h) diverges 
as t^/** at the transition temperature, in agreement with 
the exact result for the thickness of the wetting layer [T9| . 
Furthermore, the density of closed base-pairs [T^ scales as 
the surface order-parameter in wetting, {h~^) ~ \5\ [15]. 

Interestingly, it was also argued in [18] that the theory 
of critical wetting should also apply to the PBD model 
(nonzero-p). Renormalisation group analyses of three- 
dimensional critical wetting as embodied in eq. Q with 
yO = famously predict a strong non-universal critical be- 
haviour [20]. These predictions, however, are at odds with 
extensive Ising model computer simulations due to Binder 
et al. [21] as well as with experiments [22], which yield a 
mean-field-like second-order phase transition for the wet- 
ting problem. Fisher and Jin [15] suggested that this dis- 
crepancy arises from fundamental defects in the wetting 
Hamiltonian 

Hen,{p = 0) = j dx [k{Vhf + Wie-"'' + ^26-2"^] , (5) 

which should include a variable, position-dependent inter- 
facial stiffness 

k{h) = k + w'^e-"^ + w'^ahe-^'''' + ■■■ . (6) 

When supplemented with the corrected stiffness, the struc- 
ture of the wetting Hamiltonian eq. ([5]) is very similar to 
that of the PBD one eq. A more detailed compari- 

son reveals that in critical wetting the parameter w'l van- 
ishes linearly with the transition temperature and it is the 

^ Just to complete the wetting story, the Fisher-Jin improved 
Hamiltonian did not yield the desired result, namely a crossover 
to mean-field like behaviour. According to a linear renormalisation- 
group study, the presence of the term proportional to uij > is capa- 
ble of destabilising the critical wetting transition, driving the transi- 



ncxt-to-leading term, u'2e^^"'', that controls the critical 
behaviour. On the contrary, for the PBD Hamiltonian 
w'l — kp/2 > is finite and t«2 is identically zero. Thus, 
a renormalisation group analysis of the PBD model along 
the same lines as in the wetting case requires switching 
on a nonvanishing w'l and truncating the series to first 
order. Despite these differences, we shall show that one 
can avail from standard renormalisation group techniques 
developed for the wetting problem to draw conclusions on 
the order of the DNA denaturation transition. Such an 
analysis is carried out in the next section, where we prove 
that the one-dimensional melting transition for homoge- 
neous DNA sequences is continuous in (h). 

Exact decimation procedure. The wetting 
analogy allows us to perform an exact decimation 
renormalisation-group (RG) [26 analysis of the DNA de- 
naturation in the PBD model. This RG procedure was 
successfully applied to the 2D wetting transition [27] 
{Hewip — 0)) where the non-trivial fixed points can be an- 
alytically calculated for short-ranged forces [28^, and also 
for more general interactions [33]. These calculations can 
be straightforwardly extrapolated to the PBD model with 
p 7^ which, as shown before, is formally equivalent to the 
discrete version of the standard interfacial Hamiltonian for 
the 2D wetting transition because the position-dependent 
stacking interaction in the PBD model is the analogous of 
the position-dependent stiffness in the wetting case. Re- 
cent studies show that the dependence on the position of 
the effective stiffness can induce new critical phenomena, 
and it can even drive the transition first-order [30) . 

The configurational part of the PBD Hamiltonian can 
be written as H' — J^iLi where H^^^ — 

W{hi,hi+i) + {V{hi) + V{hi+i))/2. By simplicity we shall 
consider periodic boundary conditions hi = Hn+i, al- 
though this choice will not affect our conclusions. We 
define the initial transfer matrix as T'^^^hi, hi^i) = 
exp{—/3H^^^). Note that T'"^ is symmetrical under an 
exchange of its arguments. The decimation RG procedure 
defines the transfer matrix T^") at any RG step n > 1 as 

/OO POO 
dhi... d/i6_iT("-i)(&'^/i,/ii) 
-OO J —oc 

X T^''-'\hi,h2) . . .T^"-'Hh-ub'^h'), (7) 

where & > 2 is the rescaling factor and C is the wan- 
dering exponent. The normalised interaction is defined 

tion weakly first-order depending on system parameters in d = 3 |23) . 
A subsequent investigation allowed the analysis to be extended con- 
cluding that a first-order transition can appear only for dimensions 
d > 2.41 I24| . This puzzling situation was clarified only a few years 
ago by Parry et al. 1251 . who argued that the effective interfacial 
Hamiltonian for short-range critical wetting in three dimensions is 
in fact nonlocal, and that in the small gradient limit, V/i 1, it 
reduces to that proposed by Fisher et al. |23| . However, and remark- 
ably, a thorough renormalisation-group and computer simulation 
analysis of the nonlocal Hamiltonian shows no stiffness instability 
and hence the wetting transition remains continuous I25| . 
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as /3ij(") = -lnr("). Note that this functional renor- 
malisation preserves the symmetry under exchange of the 
arguments of H^"-\ This RG scheme is formally exact, 
but it cannot be solved analytically in general. Instead of 
solving numerically the RG recursive equations, we shall 
analyse the effect of the RG iterations in a subspace of the 
functional space particular, we shall consider 

the transfer matrix to be equal to the propagator corre- 
sponding to the continuum limit of the PBD Hamiltonian 
Hew = Hew[h; k, p, a, a, ^1,^2] given by eq. ^ 



T{ho, hi) = Z{ho, hi;x = 1) 



(8) 



where we integrate over all the continuum paths h(t) 
(0 < t < x) subject to the conditions h{0) = ho and 
h{x) = hi. Due to the properties of the transfer integral, 
the application of the RG scheme eq. ^ to this class of 
Hamiltonians yields 



Z'{h,h';l) = b'^Z{b^h,b'^h';b), 



(9) 



where Z'{h,h';l) is the transfer matrix associated 
with a new continuum PBD Hamiltonian iJ^^ = 
H'^^ [h'] k' , p' ,a' , a' ,w'i, W2] . The renormalized Hamilto- 
nian parameters are related to the original ones via 



k' ^ kb^^-^ , p' = p 



^b<c 



= b^c 



y- = w,.b, I = 1, 2, 



(10) 



The wandering exponent is taken as C = 1/2 in analogy 
to the wetting case, so k and p are unchanged by the 
RG iterations. On the other hand, a, a, \wi\ and W2 
increase in each RG step. In order to reveal the irrelevance 
(in the RG sense) of the nonharnionic contribution of the 
stacking interaction, we shall show that the decimation 
procedure described above can be related to an analogous 
RG scheme for a modified PB model (i.e. with p = 0). 
First, we note that the presence of the position-dependent 
term in the stacking interaction makes the definition of 
the propagator eq. ([8]) ambiguous. A similar problem is 
reported for the quantisation of classical systems with a 
position dependent mass [5^155] . We choose the following 
definition of the propagator 



Z{ho,hi,;x) = lim / dhi...dhi,^i I I K{hj,hj^i;x/b), 

(11) 



where K{h, h']x) is defined as 



K{h,h';x) 



/3fc(l -h/9e-«(''+^')) 

2'KX 



(^-^(1 + pe-"(''+'''))(/i - h')^ - xv{h, h')^ , (12) 



X exp ( -^(1 +Pe 

with a modified potential u(/i, h') given by [TOl[TT| 

v{h,h') = ^u;ie- 5 + l3w2e-''^^+^'^ 

/3W3 



+^ln(l + pe-"(''+"')), 



(13) 



where initially jSw^ = 1 . With this definition the propaga- 
tor can be understood as the result of a first RG step for 
large b and x = b before rescaling the distance (see eq. ([7])). 
On the other hand, this expression reduces to the PB case 
as p 0. The propagator Z verifies a Schrodinger-like 
equation. To obtain it, we note that for small Ax 



dZ 

Z(h, h'; X + Ax) « Z(h, h'; x) + Ax- 
ox 



dh" Z{h,h" ■x)K{h" ,h' ] Ax). 



(14) 



We expand Z{h, h"; x) and K{h", h'; Ax) for smaU \h"-h'\ 
as 

see eqs. m5\) and m6\) 

w here we define d A(/tO = ^fc(l+pe-2"'''), G{h";h',Ax) = 
v/A(/i')/27rAa;exp(-A(/i')(/i" - h')^/2Ax) and the prime 
denotes differentiation with respect to the indicated argu- 
ments. Due to the Gaussian form of G{h"; h' , Ax), we can 
evaluate trivially the integrals on h" . In the limit Ax 0, 
the resulting expression reduces to 

dZ{h,h';x) _ 1 d fldZ{h,h';x)\ 

I + w Z [W) 



dx 



2dh' V A 



dh' 



with the initial condition Z{h,h';0^) — 5{h — h'), and 
wher^ 

"•<'-)-<"''-)-^-45(fy)' 

Introducing the change of variables [34] 



77 = / dh^kih), (20) 
Z{t^,tj'-x) = A{h)-^^^A{h')-'^^Z{h,h';x), (21) 
eq. (mi) yields 

dZ{r],f]';x) 1 d'^ Z {rj , r]' ; x) 



dx 2 drj''^ 

with the effective potential v — v{r]) being defined as 



(22) 



vi-q) 



v*{h{r])) 



7 
32A 



A^ 

8A2 



h=h{rj) 



= v{h{^),h{r^))-^(^ 



(23) 



h=h(i^) 



Consequently, the propagator of the PBD model can be 
mapped onto a propagator of a PB-like model (i.e. p = 0), 
where the effect of the position-dependent stacking inter- 
action is absorbed into the definition of the variable rj and 



^ Note that this Schrodinger-like eq. (I19I I is not the same as the 
one proposed in Ref. which is in fact associated with a non- 

Hermitian Hamiltonian operator. Our approach preserves the self- 
adjointness of the corresponding Hamiltonian operator, as expected 
from the symmetric character of K(h, h'; Ax). 
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K{h'\h';Ax) « G{h";h',Ax) 



l + Axv{h\h') + ^j^^{h"~h') 



[h" - h'f I K"{h') _ 1 f k'ih') 
16 I K{h') ~ 2 \ K(hJ) 



{h" 



16Aa; 



A"{h' 



dh' 2 



A(/i') 



{h" 



(15) 



32(Aa;)2 



(A'(/i')) 



(16) 



the effective potential v. From eq. ([20)) we obtain the ex- 
pression for the variable i] 



riQi; a, k, p) = y/ /3k 



.+ i(lni±^^^^ 
a V 



2 

— 2ah 



, (24) 



so 77 ~ y/pkh for /i > a~^. On the other hand, this 
expression verifies r]{Vbh; a,k, p) = y/bri{h] \fha^ k, p). Fi- 
nally, the second term in eq. (1^51) decays exponentially at 
large distances as 



_9_ /iV 
8A V A 



2 „2„-4Qh 



2 „2„-4a/i 



8/3fc(l -Kpe-2"'') 



8^fc 



(25) 



Therefore, the effective potential v{'q) decays exponen- 
tially with 77. Substituting cq. ([2T|) into eq. ([9]), and taking 
into account eq. ([24| . we find that under a RG step the 
propagator Z renormalizes as 

Z' = Z'(i^(h-a\k',p'),T^{h'-a,k,p)-l) 

= VbZ{7]{Vbh;a,k,p),r]{Vbh';a,k,p);b) (26) 

where the renormalized effective potential parameters fol- 
low eqs. (Uni) with ( = 1/2. Note that eq. ^ is not 
the recursion relationship for 2D RG decimation scheme. 
However, as we iterate the RG equations, the variable 77 
becomes proportional to h as — >■ 0. Consequently, 
the high-tcmpcrature (HTFP) and critical (CFP) fixed 
points of the standard decimation RG procedure Tgrppp 
and T^pp, respectively, P51[?T] 



HTFP 



T, 



f^(e-^-e-^),(27) 

Bk ( fik(h-h')2 fikjh+h'f 
- e ^),(28) 



CFP 



are also fixed points for the recursion eq. (I26|) , correspond- 
ing to a^^ ~ 0. Moreover, the RG flow given by eq. (|26|) 
differs from the expression for a true decimation only in 
terms proportional to exp(— 2a/i) for large a, which lead 
to irrelevant corrections in the RG sense. Therefore, the 
RG flow close to the critical fixed point must be qualita- 
tively similar to that obtained for p — 0, and we conclude 



that the DNA denaturation transition in the PBD model 
is continuous and belongs to the 2D short-ranged, critical 
wetting universality class. 

Conclusions. — In this paper we have addressed the 
question of the order of the DNA denaturation transition 
for the PBD model. By using an exact decimation proce- 
dure, we have shown that the position-dependent stacking 
interaction is irrelevant in the RG sense, so the transition 
is continuous and in the same universality class as the 2D 
critical wetting for short-ranged forces. However, our anal- 
ysis only identifies the true asymptotic critical behaviour, 
not its range. 

If the universal critical region is narrow enough, numeri- 
cally obtained denaturation transitions in the PBD model 
may look like first-order. For typical values of the pa- 
rameters, a crossover temperature Tcross/Tm ^ 0.99 has 
been numerically estimated, Tm being the melting tem- 
perature (see Fig. 2 of Ref. 1Q\). Thus, the critical region 
turns out to be very narrow, which is a consequence of 
the entropic barrier induced by the anharmonicity in the 
stacking interaction [10 . This would explain the difficul- 
ties in determining the order of the transition that have 
been reported in the literature. 
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